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We present a study of the attractive Hubbard model based on the dynamical mean field the- 
ory (DMFT) combined with the numerical renormalization group (NRG). For this study the NRG 
method is extended to deal with self-consistent solutions of effective impurity models with super- 
conducting symmetry breaking. We give details of this extension and validate our calculations with 
DMFT results with antiferromagnetic ordering. We also present results for static and integrated 
quantities for different filling factors in the crossover from weak (BCS) to strong coupling (BEC) 
superfluidity. We study the evolution of the single-particle spectra throughout the crossover regime. 
Although the DMFT does not include the interaction of the fermions with the Goldstone mode, we 
find strong deviations from the mean-field theory in the intermediate and strong coupling (BEC) 
regimes. In particular, we show that low-energy charge fluctuations induce a transfer of spectral 
weight from the Bogoliubov quasiparticles to a higher-energy incoherent hump. 

PACS numbers: 71.10.Fd, 71.27.+a,71.30.+h,75.20.-g, 71.10.Ay 



I. INTRODUCTION 

The Hubbard model of locally interacting fermions 
plays a fundamental role in the theory of condensed mat- 
ter physics and has become a standard model to study 
correlated electronic behavior. In its repulsive version 
depending on interaction strength and doping it displays 
magnetic instabilities such as antiferromagnetism. How- 
ever, there is also evidence 1 -' 2 -' 3 -' 4 - that there is a parameter 
range where it possesses a strong instability in the pair- 
ing channel to d-wave superconductivity, which makes 
it a good candidate to describe many important aspects 
of the high temperature superconductors. Its attractive 
counterpart, the model with an onsite pairing term, has 
a simpler phase diagram, as the ground state is an s- 
wave superconductor. At half filling a degenerate charge 
ordered state can also occur. For electrons in a solid 
this model may seem inappropriate at first sight, but one 
can think of the local attraction between the electrons as 
mediated by a boson, for instance, a phonon or exciton, 
where any form of retardation is neglected.— Indeed, the 
Bardeen, Cooper, and Schrieffer- (BCS) theory for su- 
perconductivity uses a similar model with instantaneous 
local attraction albeit with an energy (Debye) cutoff. In 
ultracold atom experiments 7 - the interactions between the 
fermionic atoms in an optical trap can be tuned by a 
Feshbach resonance. For a broad resonance there exists 
a regime where the effective interaction is well described 
by a local attraction. Superfluidity has been observed in 
such systems 7 ' 8 ' 8 '' 1 ' ', also in the case where the fermions 
are confined to an optical lattice 1 - 1 -. 

When tuning the interaction in models of attractive 
fermions, such as the attractive Hubbard model, one has 
two limiting cases, that of weak coupling BCS superflu- 



idity and the strong coupling Bose Einstein condensation 
(BEC) of preformed pairs. The theoretical understanding 
which has been developed over the years is that the prop- 
erties, such as the order parameter A sc and the transition 
temperature T c to the superfluid state, are connected 
by a smooth crossover, and approximate interpolation 
schemes between these limits have been devise d 12 ' 13 ' 14 ' 15 . 
Apart from its recent experimental realization for ultra- 
cold atoms in an optical trap 7 - 8 - 9 -! -, there is experimental 
evidence that this BCS-BEC crossover has also relevance 
for strong coupling and high temperature superconduc- 
tors. It has been claimed that these superconductors 
display properties in certain parts of the phase diagram, 
such as the pseudo-gap, that can be understood in terms 
of pairs, preformed above the transition temperature T c , 
in contrast to the BCS picture, where the pairs no longer 
exist above T C M&±1 

Many aspects of the attractive Hubbard model have 
already been investigated 5 - 1 ^. However, the dynamic 
response functions have received fairly little theoretical 
attention, and it is the predictions for these quantities 
through the crossover that will be the focus of the present 
paper. One particular question concerns the fermionic 
excitations in the one particle spectral functions. These 
are dominated by sharp Bogoliubov excitations in the 
weak coupling limit. However, at strong coupling, when 
the fermions are bound to pairs, we expect a decrease 
of the spectral weight carried by the Bogoliubov quasi- 
particles. In order to investigate in detail what happens 
throughout the crossover a suitable approach to calcu- 
late dynamic quantities is required. In situations where 
the momentum dependence of the self-energy is not so 
important, such as in the Mott transition, the dynamical 
mean field theory (DMFT) has proven to be useful as lo- 



2 



cal interactions can be treated very accurately. A variety 
of methods such as perturbation theory, quantum Monte 
Carlo, as well as exact diagonalization (ED) and numer- 
ical renormalization group (NRG) are commonly used to 
solve the effective impurity model. Amongst these meth- 
ods the NRG is one of the more suitable ones to cal- 
culate low temperature spectral functions. Since it was 
originally proposed by Wilson^, it has been developed 
constantly over the years/2£ The way of calculating spec- 
tral functions has been given a solid basis by the recent 
approac h 21 ' 22 based on complete basis set proposed by 
Anders and Schiller—. So far the NRG has, however, 
not been applied to self-consistent DMFT calculations 
with superconducting symmetry breaking. Here we will 
show in detail how the method can be extended to this 
situation and present results for the spectral functions. 
Some of the results have already been published in Ref. 
I24L DMFT studies for the attractive Hubbard model 
based on other 'impurity solvers' have been carried out 
in the normal phas o 25 i 26 , and in the broken symmetry 
phas o 16 ! 27 ' 28 . There is also a recent study in two dimen- 
sions with cellular DMFT 2 ^. 

Our paper is organized as follows. The model and 
DMFT-NRG approach are described in section II. For 
this calculation the DMFT-NRG approach has to be gen- 
eralized to deal with the case of a superconducting bath. 
This generalization is described in detail in section III. 
There is a mapping from the negative U model to the 
positive one when the lattice is bipartite. In the half 
filled case this mapping can be used to check the results 
for superconductivity with earlier DMFT-NRG calcula- 
tions with antiferromagnetic order. The mapping and 
comparison of the results is given in section IV. In sec- 
tion V we compare our results for static and integrated 
quantities, such as the anomalous expectation value or 
superfluid density, with results based on other approx- 
imations. Finally in section VI we present results for 
dynamic response functions. We focus on the features in 
the one-electron spectral density. Dynamic susceptibili- 
ties calculated with the method described here have been 
reported in Ref. Hi- 



ll. MODEL AND DMFT-NRG SETUP 

The subject of this paper is a study of the attractive 
Hubbard model, which in the grand canonical formalism 
reads 



H =Y1 (**i c l°- c .?> + h - c -) ~ ^ n ™ ~ 11 



ni,in it i, 



(1) 



with the chemical potential fj,, the interaction strength 
U > and the hopping parameters Uj. c\ a creates a 

fermion at site i with spin ex, and rii^ — c| CT Ci j(7 . The 
present calculations are confined to zero temperature, 
however, an extension to finite temperature is possible. 
To study superconducting order we can include an ex- 



plicit superconducting symmetry breaking term, 

ff sc = A s ° c ^[4 )T ct_ fe4 +h.c.], (2) 

k 

with an "external field" Ag C . In the superconducting case 
in Nambu space the Green's function matrix is given by 



C -fc. I i C k.]))^ (( C -fc. I i C -M 



(3) 



where we use the notation for zero temperature retarded 
Green's functions for two operators A, B, ((A; B)) u := 

-ijdt 6»(t)e Mt ([A(t),B]) with the expectation value in 

the ground state (...). Upon including ([2j) the non- 
interacting Green's function matrix G k {u) has the form, 



(4) 



where = £k — M- F° r the interacting system we intro- 
duce the matrix self-energy S fc (tj) such that the inverse 
of the full Green's function matrix G k (uj) is given by the 
Dyson equation 



(5) 



We employ the dynamical mean field theory to analyze 
the model JTJ. As effective impurity model we consider 
the attractive Anderson impurity model in a supercon- 
ducting medium, 

#And = H^ + ^Ekcl^Ck^+^Vkiclj^+h.C.) 



(6) 



where Hi mp — J2<j £ d n o- — Un^ni with n a = d' a d a and 
d a is the fermionic operator on the impurity site. Vk 
and Afc are parameters of the medium . For the model 
([6]) the non-interacting Green's function matrix has the 
form, 



Goiuj)- 1 = oj1 2 - s d r 3 - K{uj). 



(J) 



K_{lo) is the generalized matrix hybridization for the 
medium, with diagonal part 



and offdiagonal part, 



^ 2 -(4 + a|) 



(8) 



(9) 



For a self-consistent numerical renormalization group 
(NRG) calculation of an effective impurity problem one 
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has to (i) calculate the effective impurity model param- 
eters Vk, £fc and A*, in © from a given input function 
K_{uS) and (ii) map (JH) to the so-called linear chain Hamil- 
tonian, to which the iterative diagonalization of the NRG 
can be applied. Due to the symmetry breaking the stan- 
dard formulation^ needs to be extended. The details 
of how this can be achieved are described in the next 
section. 

In the case with superconducting symmetry breaking, 
the effective Weiss field is a 2 x 2 matrix 1 (t). The 
DMFT self-consistency equation in this case reads^S 



(10) 



with fc-independent self-energy^. Hence, we use the 
NRG to solve the effective impurity problem for a given 
medium K(u>) and calculate S(w) as detailed in the ap- 
pendix A. 3. From this we can obtain the diagonal local 
lattice Green's function, which for the superconducting 
case takes the form 



p (e)(( 2 (w) + e) 



H - e][Ca(w) + e] - S 21 (w)Ei 2 (a;) ' 



(11) 

where po{s) is the density of states of the non- interacting 
fermions and = u + p — En(co>) and C2{ U} ) = u — 

P — £22 (w). The offdiagonal part is given by 



G oS {uj) 



ds-. 



p (e)E 2 i(w) 



[Cx(w) - e][C 2 (w) + e] - E 2 i(w)Ei2(w) ' 

(12) 

We denote G n = G, G 2 i = G oS and G 2 i(to) = 
G 12 {-uj)*, G 22 {uj) = -G n (-u)*. These Green's func- 
tions can be collected into the matrix G. Having cal- 
culated the local Green's function G the self-consistency 
equation (flC)|) determines the new effective Weiss field 
(w). We take the impurity model in the form (J6j) , 
and identify Gq(u>) = Q_ q (oj). Then from equation ([7]) 
we obtain an equation for the effective medium matrix 
K_{uj). In the calculations with spontaneous supercon- 
ducting order we will always consider the limit Ag C — * 
in equation J2]), where a solution with superconducting 
symmetry breaking will have bath parameters 7^ 
in the effective impurity model ((6]). In section IV we 
compare the results of our extended method with the 
ones from a well-known antiferromagnetic case in order 
to gauge the quality of the new scheme. 



III. EXTENSION OF THE NRG FORMALISM 
WITH SUPERCONDUCTING SYMMETRY 
BREAKING 

In this section we give details for the extension of the 
DMFT-NRG calculations with superconducting symme- 
try breaking. We first outline how to extract the param- 
eters of the impurity model from the medium function. 
Then we discuss the mapping to the linear chain Hamil- 
tonian with details in appendix A.l. This is a generaliza- 
tion of the scheme for the normal caseSi. In the appendix 



A. 3, we describe the generalization of the calculation of 
the self-energy via the higher order Green's functions. 



A. Parameters of the effective impurity model 

In the self-consistent procedure the parameters of the 
effective impurity model have to be determined from the 
input functions of the medium K%x and K 2 i, equations 
([8]) and |(9]). We start with the Hamiltonian in the form 
((6]) and choose a discretization in the usual logarithmic 
way to intervals 1%, J+ = (x n+ i,x n ) I~ = -(x n ,x n +i), 
x n = xoA~ n , characterized by the parameter A > 1, 
and xq large enough to cover nonzero spectral weight. 
Following the normal discretization steps^ retaining only 
the lowest Fourier component yields 

-^And = -Hirnp + tn a a,n,a a a,n,a + 7n ( a l,n,u^ 

tr,n,a a,a,n 

+h.C.) - ^n(4,™,T a i,n,| + Oa 1 «,iOa,n,t)(l3) 

We outline a procedure to obtain the parameters £™, 7™ 
and 5%. For the discretized model lfT3|) we find similar 
equations to |(8|) and ([9]), 



K ^ Z ) - 2^™ y 2 _ F a 2 ' 



K 2 i(z) = Y,<' 



Z 2 - E° 

z 2 -E« 2 ' 



(14) 



(15) 



with E% = 1/^ 2 + 6% 2 . The imaginary parts A(u) := 
— IuiKix(u) + irf)/ir and A off (w) := — lmK 2 i(u) + ir])/ir 
can be written as sums of delta functions, 

AH = ^7: 2 [<^-^n)+<* + ^)], 
n,<x 

where 

ul=l(l + B-) and vl„ = l(l i;; 



2 V E3 



2 V EB 



(16) 

with u^ a + v\ a = 1. We define the spectral weights for 
the delta function representation in the intervals I" by 

w n , a = y^ w A(w) and w n ,a = J^> A off (w). 

jo, ja 

(17) 

If we assume that E" S I" , then the equations give for 

a = +, 

= 7,t 2 "n,+ +ln X,-, (18) 
Wn,+ = Jn 2u n, + Vn,+ +Jn 2u n-Vn,-, (19) 
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and similarly for a = — . This leads to three independent 
equations to determine the four sets of independent pa- 
rameters 7^ 2 , 7,7 2 , u, h+ and u n ,-- Hence, we are free 



to choose one of them, e.g. 7„ 



+ 2 _ 



from which 



follows directly j n 2 = if n ,-. We are then left with the 
equations 



and 



«)n,+ = Wn t+ U n ,+Vn,+ + -"», 



<-), 

(20) 



(21) 



Using the equality 



(22) 



we can derive a quadratic equation for d uv>ot = u 2 a — v. 
with the solution 



2 

a, a 



*uv,-{- 



W nt+ W ni -(w nj+ W n - - 2u> 2 + ) 



+4u>n, + Wn,+\l Wn, + W n _ - W 2 _ + / 

iUn,+«V-(«V+«V- - 2w 2 + ) + wf l + + Aw^ + w n ^ + 
By definition the parameters are then obtained from 

S n — 2u na v n ^ a E n , = {u na v na )E n . (23) 

In the symmetric case, w n .+ — u>„._, this simplifies to 



l n,+ v n,+ 

such that 

W-n 




J n,+ 



■w 



2u nt+ V n + = 



(24) 



-En, 5- = 



W n .- 



-E n , G = aJl 



w 



W 



n,+ 



Apart from the condition that it lies in the inter- 
vals J", E" has not been specified, but it is reason- 
able to take a value in the middle of the intervals, i.e. 
E% = \x n + x n+ i\/2 > 0. With this choice all parame- 
ters are specified numerically and the discrete model is 
determined fully by the input functions. It can be eas- 
ily checked that this procedure simplifies to the standard 
procedure 2 ^ in the case without superconducting symme- 
try breaking. 

It is also useful to check that in the case of a mean field 
superconducto r 32 i 33 i 34 i 35 i 36 i 37 i 38 the usual expressions for 
the impurity parameters are recovered in this scheme. 
For simplicity we assume A sc <§; D in the following. Ex- 
pression (|All|) for the free impurity Green's function for 
this model yields for the medium functions analytically 
for M > A sc 



A(w) = -- 



7T yfu 



A' 2 



(25) 



and 



A ott (u;) = - 



r a s 



(26) 



With the described procedure one finds apart from a 
small correction the standard results for and 7". In 
addition we obtain 



1 + (A_H 

4 



-0(A 3 C ), (27) 

where we used an expansion both in A sc and (A — 1). 
Hence, in the continuum limit, A — ► 1, 5% = A sc comes 
out correctly as the constant mean field gap parameter. 



B. Mapping to the linear chain 

The second important step (ii) in the self-consistent 
NRG procedure is to map the discretized model (fl3|) to 
the so called linear chain model of the form, 



JV 



jV 



i?And = H imp + enfn,<rfn,(T + 2^2 ^«(/n,<r /»+!,' 



cr,n—0 
N 



+h.c.)-^A„(/l T / r t a + /„, i /„, T ), 

71=0 

with /-i,o- = da- and = >/£o, with 

n 

As usual we define the localized state 

fo,<7 = -7j= y^Xl£a+,n,<r + In 0-,fi,tr)- 

v Co „ 



(28) 



(29) 



(30) 



The orthogonal transformation between the two 
Hamiltonians needs to be more general than in the stan- 
dard case since with superconducting symmetry break- 
ing we have superpositions of particles and holes in the 
medium. We choose the following ansatz for the trans- 
formation 

fn,1 — ^ ^Q,nmQa,m.f ^Q,nm^ am ^i (31) 



and 



f n I — ^ ^o,nm^a,m,t ^ ^Q,nm^ Qjm ji (32) 



We can now derive the recursion relations for the matrix 
elements and the parameters. This is done in generaliza- 
tion of earlier work by Bulla et al.~ and the details are 
given in the appendix A.l. We find for the parameters 
of the linear chain Hamiltonian (l28l) 



£)i ^ ] (m( M a,nm v a,nm) ^^m u a,nm v a,nm, (33) 
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^■n — / $ n ( u a ,nm V a,nm) ^^ m u a,nm v a,nm, (34) 



and 



n' ,a 



Srr V ct.nn' 1 a,nn' / 1 n' \ a.rin' 1 a, nn' / 



_ F 2 fl2 _ a2 



(35) 



The recursion relations for the transformation matrix el- 
ements read 

PnU a , n +ln> = (36) 



lattice sites. The transformation from the attractive to 
the repulsive model ([38)1 yields 



c h,l 



c A,k^ + c B,fc,T' 

J _J 
c A,fc,| L B,k,T 



(40) 
(41) 



Since we assume Neel type order the quantities of the 
B-lattice are related to the A type lattice with opposite 
spin. We find 

The local lattice Green's function for the antiferromag- 
netic Green's function is obtained by fe-summation over 
the reduced Brillouin zone ^ fc — > Jde p (e)/2, 



and 



PnVa.n+ln' — (37) 



IV. COMPARISON WITH AFM DMFT-NRG 
RESULTS 

There is a canonical transformation which maps the 
attractive Hubbard model with arbitrary chemical po- 
tential to a half-filled repulsive model with a magnetic 
fields 



t 



(38) 



with q such that e tqoRi changes sign from one sublattice 
to another. At half filling the respective states with bro- 
ken symmetry, superconductivity (SC) and antiferromag- 
netic (AFM) order, correspond directly to each other. 
Hence, the quality of our new method for the supercon- 
ducting can be tested with well-known DMFT results 
from the case with antiferromagnetic orderin g 40 ' 41 . 

The mapping can be applied to map the corresponding 
effective impurity models of the two cases onto one an- 
other and we give the details in appendix B. Here we use 
the mapping (|38"1) to relate the dynamic response func- 
tions from the AFM and the SC case, and we focus on 
the integrated spectral functions for the two calculations. 
In the antiferromagnetic case in the DMFT study we usu- 
ally use the A-B sublattice basis c£ a = (c A k a , c^ B k a ), 



((cA,k,r,c\ M )) u ((cA,k,f,C i Bk1 -)) u 
((c B ,k,r> c A,k,^ u « c S,fc,T; CB,fc, T )), 



(39) 

where fe is in the reduced Brillouin zone as we have dou- 
bled the Wigner-Seitz cell in position space including two 



Ga,t,iM = - de p {e) 



(42) 



where Ca,<r(^) = OJ + p- a — Y* a , a (u>). The offdiagonal el- 
ements vanish as product of a symmetric and antisym- 
metric function, 

G A ,t,l^) = Ude po(e) j |V " 7 , o = 0. (43) 
2 J (a,i(u)Qa,i(w) - e z 



As a result, we can directly relate the diagonal local lat- 
tice Green's function Gu(oj) of the superconducting sys- 
tem to the sublattice Green's functions of the antiferro- 
magnetic system, 

G n (u>) =Ga,t,tH+ G AI,|H- ( 44 ) 
Similarly, one finds for the offdiagonal Green's function, 
Gi2(w)=G A ,t, T (w)-G A u(w). (45) 

The antiferromagnetic order parameter Aafm = Utua, 
"niA = \{ n A,] ~ n A,i), is therefore directly related to the 
superconducting order parameter A sc = [/$ , 



u 

,4> = Jduj(-hmG° s (u;; 



(46) 



The results in this section are calculated with the Gaus- 
sian density of states po(e) = e - ( £ /* ) / 'y/rrt* correspond- 
ing to an infinite dimensional hypercubic lattice. We de- 
fine an effective bandwidth W = 2D for this density of 
states via D, the point at which pa(D) = /?o(0)/e 2 , giv- 
ing D = \/2t* corresponding to the choice in reference 
I42L We take the value W — 4. 

In the following Fig. Q] we show the comparison of 
the anomalous expectation value <I> (SC case) with the 
sublattice magnetization to ,4 (AFM case). 
We can see an excellent agreement of the corresponding 
expectation values from the two different calculations in 
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FIG. 1: (Color online) Comparison of the anomalous expec- 
tation value $ in the attractive model with the local magne- 
tization m,A in the AFM DMFT calculations for half filling. 



all coupling regimes. In Fig. [2] we show the comparison 
for Green's functions for U = 1,3,6. 
We can see that for the whole frequency range the over- 
all agreement of these spectral functions is good. In the 
weak coupling case, U = 1, differences can be seen in the 
height of the quasiparticle peaks, which are sharper and 
higher in the calculation with superconducting order. In 
contrast, at strong coupling, [7 = 6, the peaks are a bit 
broader and not as high as in the antiferromagnetic solu- 
tion. It should be mentioned that for large U DMFT-ED 
calculations in the AFM state have revealed spin polaron 
fine structure in the peaks^. These have so far escaped 
NRG calculations with less resolution at higher energy, 
but improved schemes might see this in the future. Gen- 
erally, the results convey the picture of a good agreement 
for static and dynamic quantities for these two different 
DMFT-NRG calculations. 



V. RESULTS FOR STATIC AND INTEGRATED 
QUANTITIES 

Having tested the method at half filling we discuss 
results for different filling factors in this section. We 
present results for static and integrated quantities ob- 
tained with the extended DMFT-NRG method. They 
can be compared to the quantities obtained with DMFT 
calculations with other impurity solvers, like iterated per- 
turbation theory (IPT)21 or ED±£. The semielliptic den- 
sity of states with finite bandwidth 2D was used for all 
the following calculations, 
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FIG. 2: (Color online) Comparison of the spectral functions 
of SC-DMFT and AFM-DMFT calculations for U = 1,3,6 
(top, middle, bottom) for half filling. 



irD- 



(47) 



with D = 2t for the Hubbard model, t — 1 sets the 
energy scale in the following. All the results presented 
here are for T — 0. For many of the calculations we take 
the model at quarter filling (n = 1/2) as a generic case to 



analyze. For the NRG calculations we use A = 1.6 and we 
keep 1000 states at each step. In the given units U c = 2 
is the critical interaction for bound state formation in 
the two-body problem for the Bethe lattice^, and can 
be referred to as unitarity in analogy to the crossover 
terminology of the continuum system. 
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A starting point for an analysis of many quantities in 
the BCS-BEC crossover in the attractive Hubbard model 
can be mean field (MF) theory.— For a given U and fill- 
ing factor n the chemical potential /imf and the order 
parameter A sc .mf = U^mf is determined by the mean 
field eq uations. The fermion ic excitations are given by 
El = yj{e k - p) 2 + A SCjM F 2 with p = /imf + Un/2. At 
weak coupling the MF equations give the typical expo- 
nential behavior for $mf, and for large U one finds 



1. 



PMF 



-I . 



MF 



vM2 - n) 



(48) 



If p is larger than the lower band energy (in our case 
— D = —2) then the minimal excitation energy is A SCi mf 
and occurs for e k = p, which usually applies for weak 
coupling. For strong coupling and n ~ 1 the minimal 
excitation energy is also given by A SCj mf, which is of 
order U. However, for low density, n — > 0, (|48| yields 
p — > —U/2, whereas $mf and thus A SC: mf are small. 
Once p has become smaller than the lower band energy, 
the minimal excitation energy is still of order U as -E^in = 



/i 2 + A 2 C MF = U independent of n. In the low-density 

strong-coupling limit the excitation gap is given by p 
which then corresponds to the energy of the two-fermion 
bound state. 

The mean field spectral densities are given by 



MF/ 
Pk 



MF,off 
Pk 



ul5(tu-E Q k )+vl6(u + E k ), (49) 
(w) = u k v k [S(LU - E° k ) - S(cu + E° k )], (50) 



where u\ = (1 + (e k - P)/E° k )/2, v 2 k = 1 - u\. There 
are two bands of quasiparticle excitations given by ±E k , 
with weights u k for particle-like and v k for the hole-like 
excitations with infinite lifetime. 



A. Behavior of the chemical potential 

In Fig. [3] we plot our DMFT results for the chemical 
potential fi as a function of U for different densities n. 
We can see that in all cases the values tend to the mean 
field value of —U/2 for large U. The results are in agree- 
ment with the ones reported by Garg et ali^l, and as 
seen there also to the mean field values, which we did 
not include in the figure. 

In the inset we show the quantity p—Un/2, which cor- 
responds to p in the mean field theory. When the density 
is low, e.g. n = 0.15, it is seen to intersect with the lower 
band edge —2 at intermediate interactions, U ~ 3.6. 
Hence /j, plays a role to determine the fermionic exci- 
tation spectrum as discussed before. If its value does not 
change much with temperature, and fi — Un/2 remains 
smaller than —D, then no Fermi surface exists above T c , 
and the system does not possess fermionic character any- 
more as fermions are bound to composite pairs also above 
T c . For large U, \x ~ —U/2 gives the binding energy. 




FIG. 3: (Color online) The chemical potential fi as a func- 
tion of U for different filling factors n. The inset shows the 
quantity /i — Un/2. 



B. Anomalous expectation value 

One of the characteristic quantities of the supercon- 
ducting state is the presence of a finite anomalous ex- 
pectation value $. The mean field equation gives an ex- 
ponential increase for $ at weak coupling, and a quan- 
tity which only depends on the density n (|48|l in the 
strong coupling limit. In the attractive Hubbard model 
the T c increases exponentially with U and then decreases 
at strong coupling with t 2 /U due to the kinetic term 
for hopping of fermionic pairs. This is captured in the 
DMFT calculation, which investigates the transition tem- 
perature as a pairing instability from the two particle re- 
sponse function. 46 We expect the anomalous expectation 
value $ in the strong coupling limit to be reduced from 
the mean field value due strong phase fluctuations. This 
is analogous to the reduction of the antiferromagnetic or- 
der parameter in the Heisenberg model by (transverse) 
spin waves. The latter are however not captured within 
our DMFT calculations in the state with broken symme- 
try, and $ increases to a constant as in the mean field 
theory, as can be seen in Fig. [4] for quarter filling. 
The order parameter A sc dmft = C^dmft can, how- 
ever, be interpreted as a high energy scale for pair forma- 
tion then.^6 The DMFT results for <&dmft are obtained 
by integration of the offdiagonal Green's function as in 
equation (|46|) or the static expectation values calculated 
in the NRG procedure, the results of which are in very 
good agreement. MF and DMFT results show qualita- 
tively a very similar overall behavior. There is a substan- 
tial reduction of the value through the quantum fluctu- 
ations included in the DMFT-NRG result, which appear 
most pronounced in the intermediate coupling regime, 
near unitarity U c — 2. However, also at weak coupling 
there is already a correction to the mean field results. For 
instance at U = 0.7 we find $mf/$dmft ~ 2.58. This 
is comparable to the reduction found in the analysis of 
Martm-Rodero and Flores^ with second order perturba- 
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FIG. 4: (Color online) The anomalous expectation value $ as 
a function of U for n = 0.5. For comparison we have included 
the results from DMFT-IPT extracted from Ref. H3 and the 
dashed line gives the result for $mf. 



tion theory. Below U = 0.5 the ordering scale is very 
small, and we do not find a well converged DMFT so- 
lution with symmetry breaking any more. In Fig. [4] we 
have also included the results obtained by DMFT-IPT^, 
which are slightly larger but otherwise in good agreement 
with our DMFT-NRG results. 




01 2345678 



U 

FIG. 5: (Color online) Average pair density (njrij) as a func- 
tion of U for a number of different filling factors. 

for the momentum distribution is given by the integral 
over the diagonal Green's function, 

o 

n k = J da; [-ImGfe(cj)]/7r. (51) 

— oo 

In Fig. [6] we plot the momentum distribution n k calcu- 
lated from (|5Tjl in comparison with the mean field result 
for n = 0.5. 



C. Pair density 

The ground state of the system is also characterized 
by the double occupancy (n^n^) or average pair density. 
The double occupancy multiplied by U gives the expec- 
tation value of the potential energy. At weak coupling 
potential energy is gained in the symmetry broken state, 
whereas at strong coupling kinetic energy gain is usually 
responsible for Bose Einstein condensation, (n^n^) can 
be calculated directly from NRG expectation values. In 
Fig. [5] it is plotted for different filling factors for a range 
of interactions. 

In the non-interacting limit it is given by (n/2) 2 , since 
the particles are uncorrelated and the probabilities n/2 
to find a particle with spin a are just multiplied. In the 
strong coupling limit all particles are bound to pairs, and 
the pair density is given by half the filling factor, {n^n^ — 
n/2. This continuous crossover from the non- interacting 
to the strong coupling values can be seen for all densities 
with the most visible change in the intermediate coupling 
regime occurring around U c = 2. 



D. Momentum distribution 

On the mean field level the weight of the quasiparticle 
peaks is given directly by the factors u k and v 2 . as seen in 
equation (|49j) . These factors also describe the momentum 
distribution n^ F = u 2 ,. The corresponding DMFT result 




FIG. 6: (Color online) The momentum distribution calculated 
from the fe-dependent Green's function and compared with 
the MF result n™ F = v\ (dotted lines) for n = 0.5. 

For small attraction (U — 1) we can see that n k shows 
the typical form known from BCS theory dropping from 
one to zero in a small range around e k — u — Un/2. 
Therefore, some momentum states above u — Un/2 are 
occupied, but only in a small region of the size of the or- 
der parameter. When U is increased, the momentum dis- 
tribution is spread over a larger range. In the BEC limit, 
where the fermions are tightly bound and therefore very 
localized in position space, we expect the momentum dis- 
tribution to be spread due to the uncertainty principle. 
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In all cases the sum rule l/NJ2 k n k — n /2 is satisfied 
numerically within an accuracy of about 1%. There are 
visible quantitative deviation between MF and DMFT re- 
sults, but they are fairly small. Our results comparable 
well to the ones presented by Garg et al<21. 

In the experiments in ultracold gases where the BCS- 
BEC crossover is investigated the momentum distribu- 
tion can be measured quite accurately. This has been 
studied also in comparison with mean field results by 
Regal et al.— . Considering low densities for the lattice 
system, and taking into that an additional broadening 
would occur at finite temperature, a qualitative agree- 
ment of our results with the experiment can be found. 



E. Superfluid stiffness 



and the Kramers-Kronig relations for the real and imagi- 
nary parts of the Green's function such that at zero tem- 
perature D s takes the form, 

o 

D s = -^Jde k Po{e k )V{e k ) J dtu ImG£ off (w)ReG£ off (w), 

— OO 

(56) 

where G k ° (u>) is the retarded offdiagonal Green's func- 
tion (f5]). We can evaluate the expression ([56]) using the 
mean field Green's function in the form l|50p . which yields 
the somewhat simpler expression 

D 

Df F = 4 fde k p (e k )V(s k )^. (57) 



For a system in a coherent superfluid state another 
characteristic quantity is the superfluid stiffness D s . It 
is a measure of the energy required to twist the phase 
of the condensate, and therefore related to the degree 
of phase coherence of the superconducting state. Usu- 
ally, it is proportional to the superfluid density n s , which 
is experimentally accessible via the penetration length. 
Toschi et a\M have investigated the relation between T c 
and D s in the attractive Hubbard model and found that 
a linear scaling relation, as in the Uemura plot, holds at 
intermediate and strong coupling. 

D s can be calculated either from the weight of 
the delta-function in the optical conductivity or from 
the transverse part of the current-current correlation 
function^ Xjx;jx (Q, u), 



D, = D 



dia 



(52) 



The diamagnetic term Ddia is essentially given by the 
kinetic energy, 



£>dia — — -x 



d£fc po(e k )e k G k (iuj n ), 



(53) 



where G k (iuj n ) is the Matsubara Green's function. In the 
infinite dimensional limit Xjx;jx reduces to the bubble 
of normal and anomalous propagator o 16 ' 46 . From this 
and the relation -d/de k [p {e k ) V(e k )\ = Po{e k )e k and 
integration by parts one finds that the diamagnetic term 
cancels, which yields^ 



D s = 



ft 



de k po(£ k )V(e k )G% B (iLJ n )G k s (iu n ), (54) 



where V{e k ) = {At 2 - e 2 k )/3 for the Bethe lattice, 
can use the spectral representation, 



We 



(55) 



This expression can be evaluated in the limit C/ — > 0, 
A sc — > as u k vj./E k ) goes to a delta function then, and 
hence D s — > 2p (p)V(p,). 

In Fig. [7| the superfluid stiffness D s calculated from 
equation l|56p is displayed as a function of U for quarter 
filling. The dashed line shows the result as obtained from 
equation (|57|) . where the mean field Green's functions are 
used to evaluate the integrals. 




FIG. 7: (Color online) The superfluid stiffness D 3 as calcu- 
lated from the offdiagonal Green's function in equation |15G[| 
for n — 0.5. The dashed line gives the result for D a , when 
evaluated as in lf57|l . 

We can see that the results for D s of DMFT and MF 
calculation do not deviate very much. The superfluid 
stiffness is maximal in the BCS limit and decreases to 
smaller values in the BEC limit. D s is proportional to 
the inverse of the effective mass of the pairs uib ~ U /t 2 , 
and therefore expected to decrease like 1/U. The sys- 
tem in this limit consists of heavy, weakly interacting 
bosons, with less phase coherence. The results shown are 
in agreement with the ones reported by Toschi et alj^. 

Summarizing this section, we see that our DMFT-NRG 
results for chemical potential, static and integrated prop- 
erties at zero temperature are in good agreement with 
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earlier calculations based on different impurity solvers. In 
fact most of the results are in good agreement with mean 
field theory and quantitative deviations due to the fluctu- 
ations included in DMFT are not very large. One could 
therefore argue that the main features are already fairly 
well described by the simpler static mean field treatment. 
In the next section we will turn to spectral quantities. In 
contrast there certain features like the distinction of co- 
herent and incoherent excitations can only be described 
when we go beyond the mean field theory. Some of these 
extra features found in the spectral resolution are lost 
again when considering integrated quantities. 



VI. SPECTRAL FUNCTIONS 

In this section we present our DMFT-NRG results for 
the local spectral density p(ui) and the fc-resolved spec- 
tra, pk(to) — — ImGfc (u>)/7r, in the different parameter 
regimes. Before discussing these results in detail, and 
comparing them with those of Garg et al.— , we consider 
the different types of approximations used in the IPT 
and NRG calculations. These are relevant in assessing 
the two sets of results to arrive at a clearer physical in- 
terpretation. 

The approximation used in the IPT is in restricting 
the calculation to the second order diagram for the self- 
energy. This is evaluated using the Hartree-Bogoliubov 
corrected propagator for the effective impurity. If this 
propagator has a spectral density with a gap E g , then 
imaginary part of the self-energy from this second order 
scattering term can only develop in the regime |w| > 3E g . 
This threshold energy of 3E g corresponds to the min- 
imum energy for a fermion above the gap E g to emit 
a quasiparticle-quasihole pair excitation. Therefore, in 
the IPT the single-particle spectral functions Pfe(tc) have 
isolated delta-function peaks corresponding to the Bo- 
goliubov quasiparticle excitations with minimal energy 
Eg, together with incoherent continuous spectrum for 

M > ze 9 . 

In applying the NRG approach to the effective impu- 
rity problem, approximations arise in using a discrete 
spectrum for the conduction electron bath. The spectral 
functions are calculated as Lehmann sums over delta- 
function peaks, the positions of the peaks being de- 
duced from the discrete many-body energy levels and 
their weighting from the corresponding matrix elements. 
This is also the case for other methods using numerical 
diagonalization such as the ED (exact diagonalization) 
method. To obtain a continuous spectral function these 
delta-function peaks have to be broadened appropriately, 
usually with a lognormal function with parameter If 
the broadening is too large certain features blur, if it is 
too small the spectral functions has many spikes and is 
difficult to interpret. With such a broadening procedure 
it is difficult to resolve sharp features such as a gap in 
the spectrum and hence an energy E g . However, usually 
an estimate of the gap can be made when the broaden- 



ing is taken into account. For all the previous results on 
static and integrated quantities we have used a conven- 
tional broadening parameter b — 0.5, and the results for 
these quantities depend very little on the broadening. In 
this section we use smaller values in order to avoid miss- 
ing features which can be lost with the larger broadening 
parameter. 

Another aspect of the NRG calculations that can lead 
to some numerical uncertainty is in the way the self- 
energy is calculated. In equation (|A15|I it is shown how 
the self-energy can be calculated from the matrices of 
the Green's function G and the higher Green's function 
F. If one is interested in the values of lu for which the 
imaginary part of the self-energy vanishes then the whole 
expression in equation (|A15(1 has to be considered. As is 
well known for NRG calculations for the Anderson im- 
purity model^ the condition ImE(O) = for the Friedel 
sum rule can be reasonably well satisfied. However, ImE 
is never exactly zero and numerical errors can often be 
seen in small imaginary parts of the self-energy from this 
procedure. 

In Fig. [8]we present the NRG results for the local spectral 
density p(uj) and the real and imaginary parts for the 
diagonal and offdiagonal self-energies for U = 2 and U = 
4. We see that ImE and ImE are approximately zero 
for a certain range of ui for both for U = 2 and U = 4. 
ImE off is an antisymmetric function, which has peaks at 
similar position as ImE. ReE° is a symmetric function 
which for large u tends to the value A sc = £/$ of the 
interacting system l|46p and for small w can be interpreted 
as a renormalized gap. 

In the weaker coupling case U = 2 we find a trend sim- 
ilar to the IPT in that ImE deviates appreciably from 
zero only when |w| > 3E g . This can be seen on closer 
inspection of the top part of Fig. where we estimate 
Eg w 0.25, and ImE > for roughly u> > 0.75 (the wig- 
gles for smaller u> are interpreted as inaccuracies). This 
means that in the corresponding spectrum for /Ofe(w) there 
are isolated quasiparticle peaks and a continuous incoher- 
ent part in the spectrum for \lu\ > 3E g . The widths of 
the quasiparticles peaks, however, will not be precisely 
zero as in the IPT due to the very small imaginary parts. 
In the spectrum for p(oj) and U — 2 shown in Fig. [8] 
there is a sharp peak due to the quasiparticle excitations 
in p(uj) just above the gap (see also Fig. 1 in Ref. 
This quasiparticle band is very similar to the one based 
on IPT presented in Fig. 3 in Ref. [27|. In the IPT case 
there is a square root singularity at the gap edge which 
does not appear in the NRG results due to the small 
imaginary part in the self-energies. 

This picture changes in the stronger coupling case 
U = 4. Here it can be seen the imaginary parts of both 
the diagonal and off-diagonal self-energies develop a pro- 
nounced peak which falls within the region, E g < \uj\ < 
3E g . This leads to incoherent spectral weight in /Ofe(w) for 
\ui\ < 3E g . This is a difference with the IPT results where 
the imaginary parts of the self-energy are always zero for 
E g < \u>\ < 3E g and consequently there is no incoherent 
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FIG. 8: (Color online) The spectral functions, imaginary and 
real parts of the diagonal and offdiagonal self-energies plotted 
for U = 2 (top) and U = 4 (bottom) with b = 0.2, n = 0.5. 
The dotted vertical line gives the peak position of the spectral 
function, which can be roughly identified with E g . 



part of the spectrum for pfe(w) in the range |w| < 3E g . 
An explanation for this difference can be seen by exam- 
ining what happens to the local dynamical charge sus- 
ceptibility Xc{v) as U increases. Results for Xc(&) f° r 
U = 3,4,6 are shown Fig. [9l The excitation gap A c 
in this spectrum can be seen to decrease significantly as 
U increases. In Ref. [HI it was found^ that at strong 
coupling A c decreases like ~ t 2 /U. In the weak coupling 
case A c ~ 2E g , but for strong coupling E g increases with 
U, while A c decreases. Therefore, the contribution to the 
self-energy arising from the scattering with charge fluctu- 
ations can, for larger U, generate a finite imaginary part 
of the self-energy for E g < \ui\ < ZE g . The location of 
the peak in ImS appears consistent with such an inter- 
pretation. The same effect cannot arise from scattering 
with the spin fluctuations, as these have a larger char- 
acteristic energy scale (of the order of U) . For a further 
discussion of the behavior of the charge and spin gap we 



refer to Ref. 

The development of the peak in the imaginary part 
of the self-energy in the range \u/\ < 3E g leads to a dip 
in the local spectral function p(u) for U = 4 as can be 
seen in Fig. [H There is then a peak-dip-hump structure 
in p(u>) for U = 4 (see also Fig. [II]). This feature has 
also been found in calculations for the attractive contin- 
uum model^. This behavior is not visible in the IPT 
calculations (cf. Fig. 3 in Ref. H3). The most likely ex- 
planation is the restriction in the IPT to the second order 
diagram, which does not allow for any renormalization of 
the charge fluctuations. 
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FIG. 9: (Color online) The imaginary part of the local charge 
susceptibility for U = 3,4, 6 (b = 0.2). We find strong exci- 
tations corresponding to charge fluctuations there below the 
respective estimates of E g ~ 0.5, 0.9, 1.7. 

We now consider the effect of increasing U on the quasi- 
particle excitations which can be seen in the results for 
the spectral function pk{u>) for U — 1 and {7 = 4 for quar- 
ter filling shown in Fig. [lOl For both cases, {7 = 1 and 
U = 4, we can see a series of sharp quasiparticle peaks 
which are most narrow in the region Sk ~ ft, which is 
also the point where the spectral gap is minimal. In the 
latter case U = 4 in addition we find the hump with 
incoherent spectral weight as discussed earlier for p(u>). 
We have also added arrows which indicate the position 
of the quasiparticle peaks ±i?£ in mean field theory (|35)l . 
and the height gives the spectral weight. We can see that 
they describe the position of the quasiparticle excitation 
well for U = 1. For the larger U case the structure in 
the spectral function is not so well captured by the single 
quasiparticle peaks of the mean field theory. The energy 
of the quasiparticle excitations differ markedly from the 
mean field prediction and there is a significant transfer 
of spectral weight to the incoherent hump, in particular 
at high energy. The quasiparticle band (weight and band 
width) at low energy is, however, captured on a qualita- 
tive level by mean field theory. The quasiparticle peaks 
for pk {uj) always have a finite width since our self-energy 
is never strictly zero. One can infer the bands from the 
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FIG. 10: (Color online) The £fc-resolved spectral functions 
Pfe(w) for quarter filling in the BCS-limit, (7 = 1 (top), and 
towards the BEC limit, U = 4 (6 = 0.2) (bottom). The ar- 
rows show the delta-function peaks of the mean field solution 
p1(uj), where the height of the arrow indicates the weight of 
the peak. 



poles of the Green's function Ek and compare them with 
the mean field bands ±E%. For the weak coupling they 
are in good agreement. Towards the BEC limit the ef- 
fective mass tub of a boson pair is of order U. This is 
reflected in the smaller effective band width for the case 
U = 4 (Fig. QUI). The weight of the peaks in the full 
spectrum Pk(w) is in accordance with the height of the 
arrows for p k (ui). We can see that in the BCS limit (top) 
the weight in the lower band decreases rapidly to zero 
near Ek = p, whereas towards the BEC limit (bottom) it 
spreads over a much larger region which corresponds to 
what has been observed for momentum distributions in 
Fig. El 

In earlier work^l the quasiparticle properties were an- 
alyzed in an expansion around the solutions Ek of the 
equation ReGfc(w = Ek)^ 1 = 0. This led to the Lorentz- 
like quasiparticle peak of the form 



with width W(Ek) and weight w+{Ek). When a stan- 
dard broadening of b = 0.5 is used one finds a finite width 
of the peaks which increases in the crossover regime lead- 
ing to a strongly reduced life time of the quasiparticles. 
In the more careful analysis here with smaller broaden- 
ing and taking into account possible errors in determin- 
ing the self-energy we come to the conclusion that this 
is not generally correct. We should not expect a finite 
imaginary part of the self-energy to appear for u> ~ E g 
in a DMFT calculation which does not include collective 
modes^. 

Being aware of limitations in our numerical calcula- 
tions we investigate in more detail how the one-particle 
spectra near the minimal spectral gap modify from weak 
to intermediate coupling. Here we use a general scheme 
in which we analyze the peaks in the spectral function 
directly numerically and estimate the transfer of weight 
from the quasiparticle peaks to the incoherent part of 
the spectrum. We take the peak position in pk(to) for 
a given Ek as the excitation energy E^, the full width 
-Fpcak at half maximum as the width and the weight is 
determined by the integration over a region around -E£ x 
of 2F pca k. Such an analysis also applies to antisym- 
metric peak forms, and is equivalent to the other one 
for sharp Lorentz-like peaks. Note that a normalized 
Lorentz peak with width A (half width at half maxi- 
mum) integrated from — 2A to 2A yields the spectral 
weight u>2A = 2 arctan(2)/7r w 0.705. 
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FIG. 11: (Color online) The spectral functions pfe(w) for an 
£fe where the gap is minimal for quarter filling and (7=1 — 4 
(b = 0.2). The integration area, which gives the weight of the 
peaks is shown. 

We have done such an analysis for the e^-resolved spec- 
tral functions, where we consider an /,. such that the 
excitation gap is minimal. The corresponding spectral 
functions for U = 1 — 4 are displayed in Fig. QlJ We have 
included a line at half maximum for the width as well as 
marked the integration area in the low energy peak. We 
can see now very clearly how the spectrum changes from 
one coherent quasiparticle peak at weak coupling to the 
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peak-dip-hump structure at intermediate coupling. This 
is similar to what was found in the calculation for an at- 
tractive continuum models, where a sharp quasiparticle 
peak with little weight is still present at strong coupling. 
In our calculation the strong coupling limit is not easy 
to analyze as we always have some finite imaginary part 
of the self-energy leading to a finite width of the quasi- 
particle peak, which is partly spurious and tends to be 
larger at large coupling. At strong coupling the excita- 
tions occur at higher energy and we have to reduce the 
broadening further, which leads to a more spiky spec- 
trum. It should be mentioned that if the broadening is 
chosen larger (e.g. b = 0.5) then there is only one broad 
peak in the spectral function. 

The estimate of the weight of the quasiparticle peak 
Wpcak extracted by integration is plotted in Fig. [12] as 
a function of U. For weak coupling, U ~ 1, we would 
expect the mean field result = p) = 0.5. Due to 

the reduced integration range we find ui pe ak ~ 0.34, but 
division by u>2A gives a value close to 0.5 . 



0.34x- 1 1 1 1 1 

0.32 
0.3 

^0.28 \ 

(D X 

S 0.26 \ 
0.24 

0.22 X. - 

2 1 1 1 1 1 1 

1 1.5 2 2.5 3 3.5 4 
U 

FIG. 12: The weight of the peak for spectral excitation as a 
function of U for quarter filling. 



Coming from weak coupling we find a decrease as spec- 
tral weight is transferred to incoherent parts as seen be- 
fore in Fig. [TTJ This resembles the results of Ref. [U 

As discussed before due to the uncertainty about the 
imaginary part of the self-energy at low frequency the 
behavior of the width F poa k can not be reliably esti- 
mated. At weak coupling we expect the prediction of 
a real delta-function with F pca k — > to hold. Whether 
at strong coupling low enough excitations can be gener- 
ated to change this remains to be answered. We stress 
here again that the conclusions about the behavior of 
VUfmin Ek) as given in l|58|l . which was reported in Fig. 
3 of Ref. |2J, are found to be incorrect as judged by the 
present more careful interpretation of our DMFT-NRG 
calculations. 



VII. CONCLUSIONS 

In this paper we have presented an analysis of the 
ground state properties of the attractive Hubbard model 
in the symmetry broken phase in the BCS-BEC crossover. 
The emphasis has been on the evolution of spectral func- 
tions in the crossover regime. Our analysis is based on an 
extension of the DMFT-NRG method to the case with su- 
perconducting symmetry breaking. We have given many 
details of this extension in section III and the appendix. 
At half filling we have related our approach both for the 
effective impurity model and for the lattice quantities 
to earlier DMFT-NRG calculations with antiferromag- 
netic symmetry breaking. A good agreement has been 
found there, which validates the applicability of our ap- 
proach. As emphasized in Ref. H3, apart from the attrac- 
tive Hubbard model the extended method can be useful 
to study superconductivity in other models, such as the 
Hubbard-Holstein type, and also questions related to the 
microscopic description of magnetic impurities in super- 
conductors, which require self-consistent treatments. 

We have discussed our DMFT-NRG results for static 
and integrated quantities, like the anomalous expectation 
value, the double occupancy or superfluid stiffness. The 
results for these are in good agreement with earlier cal- 
culations based on different impurity solvers, and it has 
been found that most of the results are already obtained 
qualitatively well on the mean field level. 

The main interest of this paper has been to study 
the fermionic spectrum throughout the crossover regime. 
The local dynamics are very well described in our DMFT- 
NRG approach. We discussed how the behavior of the dy- 
namic self-energies changes when the interaction becomes 
larger. At weak coupling the spectrum is dominated 
by sharp symmetric Bogoliubov quasiparticle peaks as 
known from mean field theory. Contributions from 
particle-particle and particle-hole fluctuations incorpo- 
rated in the dynamic self-energies appear at higher en- 
ergy and are small, similar to those seen in the IPT ap- 
proach. However, when the local interaction is in the uni- 
tary regime and larger, the imaginary part of the dynamic 
self-energy shows a characteristic feature which gener- 
ates a peak-dip-hump structure in the spectral function. 
We argued that this features is likely to be generated by 
charge fluctuations as seen in the local dynamic charge 
susceptibility. One finds that spectral weight is trans- 
ferred into the incoherent parts (hump) of the spectrum 
on increasing the coupling. 

To answer the question whether at strong coupling the 
fermionic quasiparticles acquire a finite width one needs 
clarify over which region the imaginary parts of the self- 
energies vanish. Unfortunately, our method, in which 
spectral functions are obtained after broadening delta- 
peaks, is not accurate enough at present to allow for def- 
inite statements. It is possible that the peaks always 
remain sharp in the limit of infinite dimensions. Our 
DMFT approach does not capture spatial fluctuations 
and the gapless Goldstone mode. It would also be of 
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great interest to study how such effects give a modifica- 
tion of the discussed fermionic spectru m 48 ' 49 and possibly 
a suppression of the quasiparticle peaks. 
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APPENDIX A: NRG FORMALISM WITH 
SUPERCONDUCTING SYMMETRY BREAKING 

1. Mapping to the linear chain 

The second important step (ii) in the self-consistent 
NRG procedure is to map the discretized model lfT3|) to 
the so-called linear chain model of the general form ([28)1 , 



N 



N 



#And — e nfn t<7 fn,<T+ E Pn(fn,afn+1,. 



<r,n=0 



JV 



+h.C.) - £ A n (/+ iT /+ 4 + /n,i/n,T). (Al) 



n=0 



The orthogonal transformation has been chosen in the 
form (cf. equation ([3TJ)), 



a.,Tn 
n 



The matrix elements of the transformation obey the re- 
lations 



and 



o, 



which ensure that both operator sets satisfy canonical 
anticommutation relations. We can now derive the re- 
cursion relations for the matrix elements and the param- 
eters. This is done in analogy to earlier work by Bulla 



et ali^. We equate the representations for the media of 
(J3|) and (|Al(l and substitute the operator transforma- 
tion (|A2|) - (|A5j) . One can then read off the coefficients of 
the /^-operators (n > 0) on both sides of the equation, 
which yields 

n' ,a 
n'.a 

From this we find the expression f33|) for e n by taking 
the anticommutator with / n ,f ■ The anticommutator with 
/„ | gives expression I134II for A„. With the represen- 
tations 1|A2|1 - 1|A5|1 we can modify the equation (|A6[) to 
obtain 



(<5£, + A> 



t 

a,n',f 



ftn — 1 ^a,n- In' 

£[(A„-C) + (C + £n)v a,nn t 



Pn—X^a.n—ln' 



By comparison with lj3lj) we can read off a recursion re- 
lation for u Q ,n+in' in equation ([36| and for v a . n +i n ' as in 
equation lj37|) . The recursion relation for /3„ is obtained 
from the anticommutator of with / n +i,T which yields 

fin = ( M a,n+ln' + v a,n+ln' ) • 



With the orthonormality relations and the definitions e n 
and A„ we can find the expression in equation (|35|) . 



2. Relevant Green's functions 

In this section we briefly outline some details for the 
calculations of the relevant Green's functions and the self- 
energy for completeness^ For the Green's functions it is 
convenient to work in Nambu space, C\ — (dt,djj, with 
2x2 matrices. The relevant retarded Green's functions 
are then 



G d (u) = ((C d ;Cl), 



«d T ;<4» w ((d V ,d L )) u 



(A6) 

In the NRG approach we calculate Gn and G21 di- 
rectly and infer G22(w) = -Gu(-w)*, which fol- 
lows from G%b(u) = -Gf\{-uj) and G^ adv (u;) = 



-G 



ret / adv / 



-u>)* for fermionic operators A, B. Similarly, 



we can find Gi2(w) = G%l(— a;)*. In the derivation one 
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has to be careful and include a sign change for up down 
spin interchange in the corresponding operator combina- 
tion. 

In the non-interacting case we can deduce the d-site 
Green's function matrix of the model Hamiltonian ([6]) 
exactly. To do so we rewrite the superconducting term 
of the medium H sc by introducing the vector of operators 
and the symmetric matrix 



Cfc 



Then H sc can be written as 

H sc = ^ ' C^Ai-Ck. 



(A7) 



(A8) 



The matrix Green's function in the superconducting bath 
is then given by 5 fc («w„) = (iu n l 2 - ^fe)" 1 , 

£ fe (iw n ) _1 = iu n l 2 - e k r 3 + A fe n, (A9) 
where Tj are Pauli matrices. It follows that 

,. , iw n l 2 + e k T 3 - A sc ri 

= (i^-(4 + A|) • (A10) 

In the non-interacting case for T = 0, we have therefore 
G»- x - wl 2 - e d r 3 - £ V4 2 r 3 3 fc ^„)^3. (All) 

k 



The local full Green's function matrix G^w) -1 for the 
effective impurity model is given by the Dyson matrix 
equation 



where £(u;) is the self-energy matrix. 



(A12) 



3. Self-energy using the higher f-Green's function 

As described by Bulla et al.— there is a method to 
calculate the self-energy employing a higher .F-Green's 
function, and it can also be used for the case with super- 
conducting bath. The calculation taking into account all 
offdiagonal terms yields the following matrix equation 

G^r^M - UF{uj) = 1 2 , (A13) 

with the matrix of higher Green's functions F(u), 



- (UJ) ~ \F 21 (u) F 22 (u) 



(A14) 



We have introduced the matrix elements Fu(u>) — 
((d T nj.;d|)) w , F 12 (uj) = ((d T n x ; di)) u , F 21 (u) = 

-«d{n T ;4)>« and F 22 {u) = -{{d\ny, d;)) w . In the NRG 
we calculate Fn and F 2 \ and the others follow from 



Fi 2 (uj) = -F 21 (-ui)* and F 22 (u) = F n 
define the self-energy matrix by 

E(u;) = UFiw)^)- 1 . 



-U) 



We can 



(A15) 



The properties of the Green's function and the higher F- 
Green's function lead to the relations £12 (w) = £21 (— w)* 
and £22 (^) = — En(— oj)* for the self-energies. We 
can therefore calculate the diagonal self-energy E(w) = 
En(w) and the offdiagonal self-energy E° (u>) = £21(0;) 
and deduce the other two matrix elements from them. 
With the relation (|A15|) between G, F_ and £ the Dyson 
equation (j Al 2[) is recovered from (|A13|1 . Therefore, once 
G and F are determined from the Lehmann represen- 
tation the self-energy can be calculated from (|Al5j) and 
used in equations (HO), (HI]) and lfl2|) . 



APPENDIX B: MAPPING OF AFM AND SC 
EFFECTIVE IMPURITY MODEL 

In the DMFT calculations with antiferromagnetic or- 
dering the effective impurity model can be given in the 
following discrete form 



FM 



£n,a a a,n,a a oi,Ti,CT + Jn,a( a a,n,CT^ + h.C.) 



where we have omitted the impurity term. Notice that 
the parameters are c-dependent. In this model the sub- 
lattice magnetic order is taken to be in the z-direction, 
whereas in the model with superconducting symmetry 
breaking (JT3J) it corresponds to a transverse direction, x 
or y. Therefore we first perform a rotation in spin space 

— ( —( 



"'ix,ji,iji 4. v/2 

and also for the d-operators. This yields 



(Bl) 



HaFM = L n a a,n,a a <x,n,o + ^ ( a t,n,*da + h.C.) 

71. a, a n,a,(T 

n,a 



with 



ca A. Cot -,a 1 _,<* 

^= n,T l W , ~ o 1 (B2) 



pa _ S ",T s »,l tjaq _ '".t 



Then we do a particle hole transformation for the down 
spin similar to ([Mil. 



-4. 



(B3) 
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This gives 



One can then do a Bogoliubov transformation, 



HAFM - z2 L n{ a a,n,t a o',n,l + a-a,n,T a -a,n,i) 

- I] ^(-4,™ A + a_a,„,id T + h.c.) 



So far we have made no assumption about the parameters 
and 7" cr . In the usual scheme one has = — £" )Cr , 
such that L~" = — Z", Hence the second term in the 
first line is identical to the standard form apart from 
an additional constant, when we use the fermionic anti- 
commutation rules. In addition = £™j is normally 
satisfied, such that F" = 0. Therefore the term in the 
third line, which looks like the one for superconducting 
symmetry breaking, vanishes. We focus on the half filling 
case where one additionally has 7™* = 7"? So the other 
terms remain and one has a normal and an anomalous 
hopping term, 

H A FM = L n a Ln,v a *,n,v + ^ V n ( a a,n,<A + h.C.) 

71. a, a n,a,a 

+ W ni a i,n,A + d T a «,n,l + h.C.) 



a a,n,l 



(B4) 



to obtain the desired Hamiltonian H s A c nd in equation (f]~3)) . 
The matrix elements are determined by 



2 2 
11 — i) 



v° 



-V a W a 



v° 

* II 



^n.a u n,a. 



v« 2 + w% 2 - 

(B5) 

The parameters , 7" , (5" in (fl3l) are related to the ones 
in £Zafm by 



£ = 7n = V^+IT, (B6) 



(B7) 



We compared the numerical values obtained from the 
procedure described in section III for the SC case with 
the ones from earlier AFM calculations for half filling us- 
ing the above relations. A reasonable agreement for the 
two different calculations was found. 



1 C. J. Halboth and W. Metzner, Phys. Rev. Lett. 85, 5162 
(2000). 

2 D. Zanchi and H. J. Schulz, Phys. Rev. B 61, 13609 (2000). 

3 C. Honerkamp, M. Salmhofer, N. Furukawa, and T. M. 
Rice, Phys. Rev. B 63, 035109 (2001). 

4 M. Aichhorn, E. Arrigoni, M. Potthoff, and W. Hanke, 
Phys. Rev. B 74, 024508 (2006). 

5 R. Micnas, J. Ranninger, and S.Robaszkiewicz, Rev. Mod. 
Phys. 62, 113 (1990). 

6 J. Bardeen, L. Cooper, and J. Schrieffer, Phys. Rev. 108, 
1175 (1957). 

7 I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 
80, 885 (2008). 

8 M. Greiner, C. Regal, and D. Jin, Nature 426, 537 (2003). 

9 M. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. Rau- 
pach, A. J. Kerman, and W. Ketterle, Phys. Rev. Lett. 92, 
120403 (2004). 

10 M. Zwierlein, J. Abo-Shaeer, A. Shirotzek, C. H. Schunck, 
and W. Ketterle, Nature 435, 1047 (2005). 

11 J. K. Chin, D. E. Miller, Y. Liu, C. Stan, W. Setiawan, 
C. Sanner, K. Xu, and W. Ketterle, Nature 443, 961 
(2006). 

12 D. M. Eagles, Phys. Rev. 186, 456 (1969). 

13 A. J. Leggett, in Modern Trends in the Theory of Con- 
densed Matter, edited by A. Pekalski and R. Przystawa 
(Springer, Berlin, 1980). 

14 P. Nozieres and S. Schmitt-Rink, J. Low Temp. Phys. 59, 
195 (1985). 



15 M. Randeria, in Bose-Einstein Condensation, edited by 
A. Griffin, D. Snoke, and S. Strinagari (Cambridge Uni- 
versity Press, Cambridge, 1995). 

16 A. Toschi, M. Capone, and C. Castellani, Phys. Rev. B 72, 
235118 (2005). 

17 Q. Chen, K. Levin, and J. Stajic, J. Low Temp. Phys. 32, 
406 (2006). 

18 N. Dupuis, Phys. Rev. B 70, 134502 (2004). 

19 K. Wilson, Rev. Mod. Phys. 47, 773 (1975). 

20 R. Bulla, T. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 
395 (2008). 

21 R. Peters, T. Pruschke, and F. B. Anders, Phys. Rev. B 
74, 245114 (2006). 

22 A. Weichselbaum and J. von Delft, Phys. Rev. Lett. 99, 
076402 (2007). 

23 F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 
(2005). 

24 J. Bauer and A. C. Hewson, Europhys. Lett. 85, 27001 
(2009). 

25 M. Keller, W. Metzner, and U. Schollwock, Phys. Rev. 
Lett. 86, 4612 (2001). 

26 M. Capone, C. Castellani, and M. Grilli, Phys. Rev. Lett. 
88, 126403 (2002). 

27 A. Garg, H. R. Krishnamurthy, and M. Randeria, Phys. 
Rev. B 72, 024517 (2005). 

28 A. Toschi, P. Barone, M. Capone, and C. Castellani, New 
J. Phys. 7, 7 (2005). 

29 B. Kyung, A. Georges, and A.-M. S. Tremblay, Phys. Rev. 



17 



B 74, 024501 (2006). 

30 A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, 
Rev. Mod. Phys. 68, 13 (1996). 

31 W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 
(1989). 

32 K. Satori, H. Shiba, O. Sakai, and Y. Shimizu, J. Phys. 
Soc. Japan 61, 3239 (1992). 

33 O. Sakai, Y. Shimizu, H. Shiba, and K. Satori, J. Phys. 
Soc. Japan 62, 3181 (1993). 

34 T. Yoshioka and Y. Ohashi, J. Phys. Soc. Japan 69, 1812 
(2000). 

36 M.-S. Choi, M. Lee, K. Kang, and W. Belzig, Phys. Rev. 
B 70, 020502 (2004). 

36 A. Oguri, Y. Tanaka, and A. C. Hewson, J. Phys. Soc. 
Japan 73, 2496 (2004). 

37 J. Bauer, A. Oguri, and A. Hewson, J. Phys.: Cond. Mat. 
19, 486211 (2007). 

38 T. Hecht, A. Weichselbaum, J. von Delft, and R. Bulla, J. 
Phys.: Cond. Mat. 20, 275213 (2008). 

39 R. Bulla, T. Pruschke, and A. C. Hewson, J. Phys.: Cond. 
Mat. 9, 10463 (1997). 

40 R. Zitzler, T. Pruschke, and R. Bulla, Eur. Phys. J. B 27, 
473 (2002). 

41 J. Bauer and A. C. Hewson, Eur. Phys. J. B 57, 235 (2007). 

42 R. Bulla, Phys. Rev. Lett. 83, 136 (1999). 



43 G. Sangiovanni, A. Toschi, E. Koch, K. Held, M. Capone, 
C. Castellani, O. Gunnarsson, S.-K. Mo, J. W. Allen, H.-D. 
Kim, et al., Phys. Rev. B 73, 205121 (2006). 

44 A. Martin-Rodero and F. Flores, Phys. Rev. B 45, 13008 
(1992). 

46 C. A. Regal, M. Greiner, S. Giorgini, M. Holland, and D. S. 
Jin, Phys. Rev. Lett. 95, 250404 (2005). 

46 T. Pruschke, D. L. Cox, and M. Jarrell, Phys. Rev. B 47, 
3553 (1993). 

47 P. Pieri, L. Pisani, and G. C. Strinati, Phys. Rev. B 70, 
094508 (2004). 

48 K. Borejsza and N. Dupuis, Europhys. Lett. 63, 722 

(2003) . 

49 K. Borejsza and N. Dupuis, Phys. Rev. B 69, 085119 

(2004) . 

60 R. Bulla, A. C. Hewson, and T. Pruschke, J. Phys.: Cond. 
Mat. 10, 8365 (1998). 

61 These results are robust with respect to broadening as the 
excitation can be seen directly in the raw data. 

62 A feature of the infinite dimensional model is that it does 
not include a collective Goldstone mode. A coupling of the 
fermions to the Goldstone mode in a more general model 
can lead to a damping of low energy quasiparticles 



